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Abstract 

This paper proposes a solution to Stokes' paradox for asymptotically 
uniform viscous flow around a cylinder. The existence of a global stream 
function satisfying a perturbative form of the two-dimensional Navier-Stokes 
equations for low Reynolds number is established. This stream function 
satisfies the appropriate boundary conditions on both the cylinder and at 
infinity, but nevertheless agrees with Stokes' original results at finite ra- 
dius as the Reynolds number tends to zero. The Navier-Stokes equations 
are satisfied to a power-log power of the Reynolds number. The drag on 
the cylinder is calculated from first principles and the free parameter of 
the approach can be chosen to give good agreement with data on drag. In 
this revised working paper we put our approach on a firmer mathematical 
basis using the Helmholtz-Laplace equation as a linear approximation to 
the Navier-Stokes system. In so doing we demonstrate the instability of 
the original paradox. We also demonstrate the absence of a paradox of 
Stokes- Whitehead class, and give further theoretical constraints on the 
free parameters of the model. 

Key Words: Stokes Paradox, Fluid dynamics, Stokes flow, Stream function, 
Biharmonic equation, Helmholtz equation, Low Reynolds number 



1 Introduction 

The difficulty in establishing a sensible global solution to the problem of low R e 
(Reynolds number) viscous flow around simple objects, where the flow is uniform 
at infinity, has fascinated applied mathematicians for just over 150 years. Stokes 
(1851) established that there was no solution to the two-dimensional, steady, 
incompressible, Navier-Stokes equations for asymptotically uniform flow around 
a cylinder in the biharmonic limit. This situation is now routinely described in 
modern tutorial discussions. See, for example, Chapter 7 of Acheson (1990) for 
an exercise on Stokes' paradox and a discussion of the corresponding situation 
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for the sphere. The biharmonic equation can be solved in a neighbourhood of 
the cylinder (i.e. a circle) by any stream function of the form 

ip = Csin 8(r log r — —r + —) (1) 

but there is no choice of C for which if> ~ rsin(9 for large r. The analysis of 
this problem has lead to several classic papers (Oseen, 1910; Lamb, 1911) and 
its understanding through the use of matched asymptotic expansions (MAE) 
is one of the triumphs of perturbation theory. The reader is referred to Van 
Dyke (1964) for his classic survey of the work of Kaplun (1957), Proudman and 
Pearson (1957) and other key references. 

However, the MAE approach, despite its immense power and diversity of 
expanding applications, does not give a clean resolution of the original difficulty 
in that such methods rely essentially on computing and then matching solutions 
to the problem defined on two regions: close to the cylinder and far from the 
cylinder. The purpose of this paper is to address the problem of finding global 
solutions for the low Reynolds number limit, i.e. to resolve the original paradox. 
There are, of course, other approaches to the paradox. Recently, Villas Boas 
|llj has considered the problem from a three-dimensional perspective and points 
out that there is then no paradox. 



2 Viscous incompressible flow in 2D 

A large class of fluids can be characterized by their density, p, a scalar field not 
presumed to be constant, and their dynamic viscosity p. The flow is charac- 
terized by a velocity vector field v, and an associated scalar pressure field p. 
Conservation of mass is expressed by the continuity equation 

^ + V.(pv) = (2) 

and the conservation of momentum is expressed by the Navier-Stokes equation^ 

di) 

p(^+w-Vi)) = -Vp+/iV 2 ji (3) 

If the fluid is incompressible in the sense that p is a constant in both time and 
space, we have the condition: 

V.v = (4) 
To analyze matters further, we introduce the vorticity vector 

w = V x v (5) 



1 Here V 2 acting on vectors should be understood as the ordinary Laplacian acting on 
Cartesian components. 
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In the following discussion wc demand incomprcssibility but allow for non-zero 
vorticity. Using simple identities from vector calculus the Navier-Stokes equa- 
tions may then be recast in the form 

du 1 

P (^-VXLl) + V(p+ -pv 2 ) = -LlV XLJ. (6) 

Taking the curl of this, we arrive at the vorticity equation 

^+«-Vw-w.V« = i/V 2 w (7) 
at 

where the kinematic viscosity v — \ij p. 



2.1 The stream function 

Since the velocity field is divergence-free, we may introduce a vector potential 
such that 

v = V x vj/ (8) 
and furthermore we may choose it so that it is divergence free: 

V.* = (9) 

The vector potential can be reduced to a single function when there is an ap- 
propriate symmetry. The resulting object is a stream function. For example, 
planar 2D flow is obtained by setting (and note that this automatically satisfies 
satisfies the divergence condition) 

-^{x,y,t)e z (10) 
Next we note that under the assumption that satisfies V .}&_ = 

u= -V 2 * = e z V 2 <f (11) 
and the vorticity equation becomes 

For problems where it is possible to identify a natural length scale L and 
a natural speed U, it is normal practice to perform a non-dimensionalization 
of the variables and introduce the Reynolds number R e = UL/u. Then the 
Navier-Stokes equation becomes (after rescaling the variables suitably to ip, x, y, r): 

VV-*e|-VV = i?e*^ (13) 

ot o{x,y) 

Throughout this paper we shall work in units in which the radius a of the 
cylinder is taken to be unity. In such units the Reynolds number is here based 



3 



on the radiu^jand is given by R e = U jv. The two stream functions are related 
by if? — Utp and furthermore X = x, Y = y etc. 

The old historical approach to the limiting case when R e — ► and the flow 
is time-independent is to take the view that the non-linearities may then be 
ignored (provided the non-linear term is well behaved) and the time-independent 
Navier-Stokes equations reduce to 

VV = (14) 

which is the biharmonic limit, also known as Stokes flow. It is now well known 
(see for example, Chapter 8 of Van Dyke (1964)) that the neglect of the non- 
linear terms can lead to inconsistencies, as is evidenced by the lack of any 
solution for asymptotically uniform two-dimensional flow past a cylinder. 

Here we introduce the scalar vorticity function u> = V 2 V> and write the 
Navier-Stokes as the pair 

uo = v 2 v> 

du d{M (15) 

V OJ - H e — — H e — 

dr d{x, y) 

If we had exponential growth in vorticity, lj — e^ T Cj(x,y), then the linearized 
form would be 

2 1 (16) 

V 2 UJ - R e fIUJ = 

or in terms of a single condition: 

VV - R e ^ 2 ^ = , (17) 
which is the Laplace-Hclmholtz equation. 



2.2 The instability of Stokes' paradox 

Let us focus temporarily on the linearized time-dependent case based on the 
Laplace-Helmholtz equation. We let e 2 = i? e /i. It is an elementary exercise 
(and we shall give equivalent details in a different context later) to establish 
that for all e > 0, the Laplace-Helmholtz equation 



vV - e 2 v 2 v> = o 



(18) 



has a solution of the form 



r- 1 



(e) 
eK (e) 



2K X (re) 
sK (e) 



sm(6) 



(19) 



2 it is also common to base R e on the diameter. 
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that satisfies 



v dr (20) 

i/j ~ e pr rsin(0) as r — ► oo 

In other words, no matter how small the value of /i, the paradox does not exist. 

To summarize: If we have some time-dependence in the simplified form 
of exponential growth in vorticity the linearization is precisely the Laplace- 
Hclmholtz problem rather than the biharmonic problem. Addition of this type 
of extra variability, no matter how small, shows that the paradox evaporates, i.e. 
the existence of the paradox is unstable. Villas Boas [TT], who has considered 
the problem from a three-dimensional perspective, also demonstrates that the 
paradox evaporates when extra spatial variability is incorporated. 

In subsequent analysis we will use the Laplacc-Hclmholtz model as a ap- 
proximation for the time-independent 2D case, but where the modification is 
regarded as a simple way as approximating the combined effect of the non- 
linearities via a linear term. The question is as to whether we can do this in a 
sensible and self-consistent manner. We shall look at this in a variety of ways, 
starting with a rather ad hoc approach. 



3 Reorganizing the low Reynolds number 
Navier— Stokes equations 

The production of a global solution to the Navier-Stokes equations requires a 
slightly unusual approach. We take the time-independent equation 

as our starting point and note some obvious facts. First, if we have a ipo satis- 



fying Laplace's equation, then it satisfies equation Eq. (21 1 identically. Second, 
in the high Reynolds number limit, although this is a singular limit, we note 
that the Jacobian term should then vanish identically, so that any solution of 

VV = fW) (22) 



satisfies Eq. (21 1 when R e — > oo. This will be true of course in the linear case 
when f(ip) oc ip and this last equation is of Hclmholtz type. Now let A be a 
(possibly complex) number of order 1 and /? a parameter to be determined. In 
particular, if we consider solutions, ipi, of a parametrized Hclmholtz equation 
in the form 

V 2 Vi = A.Rf (23) 
then the Jacobian will vanish and we note that 

vvi-^ yy ^^i (24) 

d{x,y) 
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Next suppose that we consider a more general stream function ip as 



i/> = V>o + "01 (25) 
Then a trivial calculation tells us that 

vV - i? e ^? = x 2 Ri ipi - xRi^ d -^tA (26) 

9(z,y) <9(a;,y) 

Our approach is therefore to combine solutions for potential flow, i/jq, with a 
solution ipx of the Hclmholtz equation Eq. (23) that also satisfies the fourth 
order equation given by Eq. (24) without the identically zero non-linear term: 

VVi = A 2 i?f Vi (27) 
When i? e — > this PDE of course approaches the biharmonic equation provided 



(3 > 0. Then Eq. (26 1 offers the possibility that we can solve a PDE of the form: 



?4( D d(v>,v 2 v) 



W - fle g ' " ; = 0{Rt)f{x,y) (28) 
where fc and / can be calculated and analyzed. The introduction of the term 



of order R^P in the right side of Eq. (26 1 may seem like an artificial device, but 



given that (a) we are seeking solutions of the problem as R e — > 0; (b) we are 
not modifying higher derivatives in the Navier-Stokes equations; (c) we have a 
valid perturbation equation in Eq. (28), we consider that it is valid to proceed 
with this approach. 

3.1 Other justifications 

The author recognizes that the argument in the previous sub-section is some- 
what ad hoc, and initially gives us no idea how to justify the choice of A or (3, 
or, as we shall introduce, the relevant composite parameter e. However, there 
is another other route to justifying this approach. Let to = V 2 V>- Then the full 
2D time-independent incompressible Navier-Stokes equation is precisely 

V 2 u;-i* e f^ = 0. (29) 
d{x,y) 

We already know from the existence of the Stokes paradox that simple lineariza- 
tion by setting R e — is hopeless. We might therefore consider starting from 
other linearizations, for example, the Helmholtz linearization 

V 2 lu - e 2 uj = , (30) 

where e has to be determined cither from experiment or from deeper theoretical 
considerations. The linear starting point for analysis is then not the biharmonic 
system but the Helmholtz-Laplace equation: 

V 4 V - e 2 VV = . (31) 
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In principle then we can consider justifying the choice of e by matching e 2 uj 
with R e ^ (t'y) as closely as possible under a suitable norm. We will prove in the 
next section that for all e > the Helmholtz- Laplace equation admits a solution 
satisfying both the boundary conditions on the circle and the velocity condition 
at infinity, obtained by considering elementary solutions that are linear com- 
binations of solutions of the Laplace equation and solutions of the Helmholtz 
equation. 

Note that this approach and our ad hoc scheme are linked by the relationship 

e 2 = \Rf (32) 

The determination of A and (3 are interesting challenges. We shall sec in Section 
8 that there are good theoretical grounds for setting (3—1. For now we will 
leave both parameters general. 



4 The case of asymptotically uniform flow past 
a circle 

We now consider the well-trodden route to the analysis of a stream function asso- 
ciated with a flow that is uniform at infinity and satisfies appropriate boundary 
conditions on r = 1. Using polar coordinates, we therefore want 

tp(r, 9) <~ rs'md as r — > oo (33) 

and 

^, = = ^ on r = 1 (34) 
or 

We build the solution for ipo and "01 a $ follows, under the assumption that ipo 
satisfies the Laplace equation and tpi the Helmholtz conditions. Any sum of 
the two will satisfy the Laplace-Helmholtz equation. Given that the angular 
behaviour at infinity is oc sin 9 we make the standard assumption and seek 
solutions 

4>i = 9i ( r ) sin 9 for i = 0, 1 (35) 
The potential flow part (as usual) will be taken to be 

5o = 1 + B/r (36) 

for some constant B. Our analysis will differ from Stokes' classic (Stokes, 1851) 
treatment in that gi does not satisfy the separated biharmonic equation. Instead 
we use the Helmholtz equation, which in separated form is just 

9» + ViW-(e 2 +^ 5l (r)=0 (37) 

where e = VXR^. The general solution to this ODE is given by 

ffi (r) = aK^er) + (3h{er) (38) 
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where h,Ki are standard "modified" Bessel functions. It is now quite clear 
that we can construct g\ so as the preserve the boundary condition at infinity 
by setting j3 = 0. The K\ function decays exponentially if e > 0. We are left 
with two arbitrary constants that can be determined by satisfying the boundary 



conditions on the circle given by Eq. (34 1. This is a matter of elementary algebra 



using some standard Bessel function identities. The final result for the total 
stream function can be simplified to: 



4, 



l 



2gi(e)\l 2K 1 {rt) 



eK (e)Jr eK (t 



sin(0) (39) 



We note some interesting facts about this expression. First, it has the right 
behaviour as r — > oo provided e > 0, as the Bessel function of r decays expo- 
nentially Second, if we keep r fixed and finite and let e — ► + we obtain 

sin 6 ( , 1 1 \ 

V ~ \ T7TT~\ \ T log r - ^ r + (40) 

log(2/e)- 7 \ S 2 2rJ 1 1 

and we recover a multiple of Stokes' (1851) solution satisfying the boundary con- 
ditions on the sphere (but not at infinity), and, furthermore, the multiple is now 
reminiscent of that arising from the methods of matched asymptotic expansions. 
The claim therefore is that it is Eq. ( |39| that essentially resolves the paradox, 
as all the boundary conditions are satisfied, but the low Reynolds number limit 
for finite r does not. The third observation is that in the neighbourhood of the 
surface of the sphere, 

^sin^r-l) 2 ^ (41) 
and that for small e this is given by 

(r-l) 2 ain(g) ( 2) 

-log( e ) + log(2)- 7 + °^ (42) 

We also note that the limit as e — > oo is just the potential flow limit: 

V>~(r--)sin0 (43) 
r 

Mathematica code for the evaluation of the stream function and velocity field is 
given in the Appendix. 



4.1 Choice of A and (3 

At this stage we have no knowledge about how to fix the parameter e = VAiif , 
and we now turn our attention to this issue and a more detailed analysis of 
Eq. ( 26 1 . We note first that our working assumption of a power dependence of 
e on R e has not actually been necessary. All we need is that e — * 0+ as R e — > 0. 
We also note that the choice (3—1 and A = l/(4e) recovers the result (Van 
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Dyke, 1964) that the limiting stream function for fixed r and small R e is, from 
Eq. gO I 

sin ^ ( i 1 1 A 

w ~ ; — ; ; i r log r r H (44) 

^ log(4/i? e ) + i - 7 V ^ 2 2r^ 1 ; 

so that the low R e limit matches exactly the first term of the MAE result]^] 
Clearly other choices of A can be considered, as can other powers of R e or even 
a more general e still. However, in Section 8 we will establish a theoretical basis 
for estimating these parameters, and argue that (3 = 1. 



5 Analysis of the remainder 



We must now analyze the right side of Eq. ( 26 1 we call it T and we do so 



without any assumption as the the form of e. We have, changing variables in 
the Jacobian to polar coordinates: 



eV - Re 



A r 



dipo d^i dtp! dipo 



dr 89 



dr 86 



(45) 



Recall that we set ipi = gi(r) sin 6. Having imposed the boundary conditions, 
we have 

_ 1 / , 2K x {e) \ _ _ 2K x {er) 



9o 



9i 



It follows that 

T = e 4 smdgi - R e e 2 sin 6 cos 6 J(e, r) = T HA +T NL 
where J is the reduced Jacobian: 



1 f dg dgi 

J = - \ 9i^ 9o-^- 

r V or or 



(46) 



(47) 



(48) 



and Tha — e 4 sm8gi denotes the 'Hclmholtz artifact' introduced by our ap- 
proach, and Tjvl denotes the non- linear Jacobian term. Some calculation with 
Bessel function identities leads to 



J(e,r) 



2(K (e)K 2 (re)-±K Q (re)K 2 (e)) 



(49) 



from which it is manifest that J(e, 1) = and the non-linear terms vanish on the 
cylinder boundary irrespective of the choice of e. Furthermore, the asymptotic 
behaviour of the Bessel functions tells us that provided e > 0, 



J(e, r) — 0(^=exp(—er)) 



(50) 



3 for analyses basing R e on the diameter we have log(8/i? e ) etc. in the denominator. 
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as r — > oo. We deduce that Tnl is a bounded function for all r. What is its 
order as a function of R e l A straightforward estimate may be given by looking 
in a neighbourhood of the cylinder. It is straightforward to establish that 

dJ ( 2K 1 (e) \ 2 4 

dr^ 1 \K (e)J ~ e 2 (-log(6) + log(2)- 7 )2 1 > 

where the latter expression applies as e — > 0. We deduce that in the immediate 
neighbourhood of the cylinder 

r - = °(o^) < 52 > 

which is o(R e ) provided only that e — > as J? e — > 0. Although J grows to 
a maximum as r is increased away from unity, before reaching a maximum 
and then decaying for large r, some numerical experiments with Mathematica 
confirm that the maximum of e 2 J slowly decreases as e decreases to zero. So 
we can assert that Tnl is well behaved and is o(R e ). When we consider Tha, 
it is easy to see that g\ has a maximum on r = 1, and that 

T H A\ r =i = e 4 sme 9l (l) = 0(e 2 /loge) (53) 

as e — > 0. So this term also behaves. If we desire that the Helmholtz artifact 
tends to zero faster than the non-linear term (which is desirable for the credi- 
bility of our approach) it is then natural to specify e in the form V\R^ and to 
demand that f3 > 1. In particular, the choice (3 — 1 arranges that 

T HA = 0(R 2 J log R e ) , T NL = 0{RJ log 2 R e ) (54) 

and we have established that T — o(R e ). It also decays exponentially at infinity. 



6 The drag on the cylinder 

The calculation of the drag can be done by integrating the pointwise force on the 
cylinder over its surface. The pointwise force has two components. One involves 
the local rate of strain in the fluid, and the other is the pressure force. The first 
requires a purely local calculation, but the second requires an integration of the 
pressure equation from infinity to the cylinder. A question is how this second 
part can be carried out without any global representation of the flow field. For 
the specific case of low Reynolds number calculations with a certain symmetry 
there are ways around the problem that we shall exploit presently. To make 
these matters clear, we shall summarize a first-principles calculation of the force 
using the results for cylindrical polar coordinates for the rate of strain as given 
by Appendix A of Acheson (1990). The fluid velocity is given (in our units) by 

^^w^^ (55) 
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On the cylinder the rate-of-strain tensor has components e rr — = egg, and 

r d 2 ip 



2e r g = U 



dr 2 



(56) 



The stress tensor = — pSij + 2paj , so the force on the cylinder boundary is 



t = T.e„ 



J. TT^ 



(57) 



We project this into the x-direction and integrate over the circle to get the 
following formula for the drag D (force per unit length on the cylinder) : 



D = 



i-2-rc 

Jo 



d 2 ib 

p cos 8 + aU—^r sin( 
or 1 



d6 



(58) 



For low R e we estimate the pressure using the momentum equation in the form 

Vp = -fj,V A (V A u) (59) 
and a short calculation gives 



Vp = pU 



-^(VV) & + |;(VV) a 



(60) 



We need to integrate this in from infinity to a general point on the circle r — 1. 
In general this is awkward without a global form of ip. But if ip has the form 
V> = g(r) sin 6 then we have 



19 / 

V 2 ?/> = sin6*£[g(r)] where C[g(r)] = -— ( r 
and the pressure equation becomes: 



dg 
dr 



.9 



Vp 



^ cos 6>£[ 5 (r)] + e e sin 0— £[ 3 (r) 



(61) 



(62) 



We do the integration from infinity along 8 = n/2 and then work around the 
circle r — 1. Carrying this out we obtain 



d 

p(l,8) =p 00 - fj,U cos 8— £[g{r)}\ r=1 
We now evaluate the total integral for the force per unit length to obtain 



D = jiU-K 



(63) 



(64) 



The drag coefficient for such a 2D problem is defined following Tritton (1959) 
as Cd = \Do/(pU 2 )\ and is now easily seen to be given by 



C D = 



R„ 



(65) 
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evaluated on r = 1. Given the vanishing of g(l) and g'(l) this simplifies further 
to 

^ = ^|.g (3) (l)l (66) 

In the case given by the MAE approach the stream function is given in the 
neighbourhood of the cylinder by Eq. (1) and then g and Cd are given by 

9 (r) = c(rlo g r- r 2 + l), C D = ^ (67) 



In our new model g as is given by the radial part of Eq. ( 39 1 , and this time, 
using some Bessel identities, we obtain 

We recall now the original assumption that e = \/Ai?f . If we fix /3 = 1 but 
allow A to vary the drag coefficient of our model is then 

C D = 2nXR e K ^ Re) (69) 

We can plot the drag coefficient from our model, assuming that (3—1, with 
various choices for A, and compare the results with those of Tritton (1959), 
which used an R e based on the diameter. Tritton 's data has been converted by 
halving his Reynolds number. Also, in making a comparison with Fig. 8.5 of 
Van Dyke (1964), it appears that Van Dyke has plotted Cd/{^k), as otherwise 
it is not possible to reconcile that plot with Tritton's data. In Figure 1 we show 
the data for the first three fibres used by Tritton (those that consider the lowest 
R e ) and 

• our model with A picked to match the MAE result, A = l/(4e) — 0.09197; 

• our model with a least squares best fit A ~ 0.04452; 

• the one term MAE result based on C = Ai(R e ) = (log(4/i? e ) -7+1/2) -1 ; 

• the two term MAE result based on C = Ax(R e ) - 0.8669Af (R e ). 

The choice of constants in the MAE approach is itself somewhat arbitrary. That 
made by Kaplun (1957), giving Eq. (40), is little more than convenience. These 
results give support to the new model and the agreement with experiment we 
get by taking A about one half that implied by the MAE approach is rather 
tantalizing. 
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Figure 1: Comparison with Tritton (1959) data 

7 Iteration: a Stokes- Whitehead anomaly? 

Another use for the type of low Reynolds numbers solution developed here is 
to provide a basis for iteration of the solution in powers of R e . The presence 
of the paradox obstructing a global representation in the standard approach 
makes this impossible. Furthermore, we know from a corresponding analysis of 
the 3D spherical problem that even when the base solution makes sense, even 
the first iteration may fail: the Stokes- Whitehead paradox emerges. So it is of 
considerable interest to investigate the iteration of the system. In doing so we 
will keep e general initially. 

We summarize the solution developed thus far. With tp as given we have, 
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combining equations (26 45) and our definitions , 



v^-R e ^^ = T = r HA + r NL , (70) 

d(x,y) 



where 



and 



2Ki (er) 

T HA = e 4 sin 9l = e 4 sing U / = lHA (e, r) sin 9 (71) 



Tml = — R e e 2 sing cos 9J(e,r) 

p 2 • fl fl 2(K (e)if 2 (re)-^if (re)if 2 (e)) 

= -fi e e smflcosfl— — ^ (72) 

^o(e) 2 

= 7jvl(e 5 r) sin cos 

where these equations also serve to define the radial terms ^xy- From now on 
we will work in terms of these radial functions, so it is useful just to write down 
their explicit forms: 

7M ( ei r),2 £ 3 |M (73 ) 



7A r -t( e ! r ) = ~2i? e e if (e) 2 



7.1 Development of the iteration 

We now attempt to refine the solution by writing the total solution as 

* = V + V-2 (75) 
Ideally we would like to arrange that 

o = v %- fl .ffi^B (76) 

5(x,y) 

The right side of this is given by the expansion 

4 ^,v 2 ^) 4 a^vV) ^,v 2 ^ 2 ) vy 2 ) 

V lp — H e — ; r h V tp2 — ti e 7T, s tC e — ; r H e — r 

d{x,y) o{x,y) o{x,y) o{x,y) 

(77) 

and we write the remaining parts of the Navier-Stokes equation as 

=V 4 i/> 2 + 1nl(c, r)sm8cos6 + jHA(t,r) sin# 

g^VV) R d{^^ 2 ) R d{^ 2l V 2 ^ 2 ) (78) 
d(x,y) e d(x,y) e d(x,y) 

We do not know how to solve this full non-linear system, and propose instead 
to treat a linearized form. The question now is whether to work with the raw 
system 

= V 4 V> 2 + n r) sin 9 cos 9 + jHA(e,r) sin (79) 
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or to propose a Hclmholtz-dampcd variation along the lines of our original 
approach. We also need to understand which non-linear term to treat first. 
Provided e is of the form e = V\R e , for example, we know that Tjvl is the 
lowest order in R e : T NL = 0(R 2 e / log 2 R e ), so we will give priority to killing 
this lowest order correction. 



7.2 Naive analysis of the raw forms 

To develop the solution we need to identify the Greens' s function for the bihar- 
monic operator. In fact we do not need the full form as the angular structure 
of the right side is simple - we need only look for a pair of appropriate radial 
Green's functions satisfying appropriate boundary conditions on the cylinder 
and at infinity. To this end we write 

■0 2 (r,9) = V3 (r) sin 9 cos 9 + V>4 (r) sin 9 (80) 

and require that 

V 4 (0 3 (r)sin6»cos6>) = -^ NL (e, r) sin 6 cos 6 (81) 

and 

V 4 (04 (r) sin 0) = -7h A{e,r) sin9 (82) 

We also require that tpi vanishes as r — > 00 and both ipi and its first derivative 
vanish on r = 1. To proceed further we recall the elementary form of the 
Laplacian in 2D cylindrical polar coordinates: 



So if / is of the form / = h(r) sin(fc#) then 



Define the fc'the radial modal Laplacian operator as 

k dr 2 r dr r 2 
We now have a pair of ODEs in the form 

*%Mr) = -lN L (e,r) (86) 
C\Mr) = ~lHA{e,r) (87) 

In a simpler problem we would now construct a pair of Green's functions Gi(r, x), 
i = 1,2, with the properties that 

£ 2 i G i {r,x) = 5{r-x) (88) 
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subject to, if it is possible, 



lim G t (r,x) = (89) 

r— >oo 

Gi(l,x) = (90) 

^Gi(r,x)\ r=1 =0 (91) 

However, there are significant complications. The Hclmholtz term we have 
thought of as an modelling artifact has first to be understood as an optimal 
linear approximation to the full non-linear system. The first steps on this are 
considered in Section 8. We will eventually have to manage the fact that no 
G\ function exists satisfying these boundary conditions. However, it is now 
consistent to proceed to construct G 2 , as this (a) exists, and (b) is needed for 
the management of the full non-linear system and its solution. 

7.3 Integration of the inertia term 

We will now focus on the lowest order correction arising from the inertia term, 
for which G 2 is the required radial Green's function. The building blocks for 
this are two different solutions of the ODE 

d 2 Id 2 2 \ 2 

.is + ?!4)«"- - (92) 

The solution to the homogeneous problem is of the form 

<f>(r) = ar 4 + br 2 + c + (93) 

To try to build a Green's function, we would write down a pair of solutions: 
<t>L{r) = a L r 4 + b L r 2 + cl + ^|, 1 <r < x , 
<t>R.(r) = a R r A + b R r 2 + Cfl+^f, x<r<oo. 

The condition at infinity requires that we set a R = = b R , so that 

dfi 

4>r{t) = c R + x < r < 00 . (95) 

The boundary conditions on the cylinder require that 

= a L + b L + c L + d L 
= Aa L + 2b L - 2d L 

So the inner function is given by 



(94) 



(96) 



Mr) ^ t ^ + a L r 4 + b L r 2 - 2b L - 3a £ (97) 
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We now have four unknowns that must be obtained by imposing junction con- 
ditions to obtain the delta-function. The junction conditions are: 



4>r{x) = <Pl{x) 
<t>' R {x) = <t>' L {x) 

4>'k{x) = 4>l{x) (98) 

<t>f{x) = <t>f{x) + l 

These equations were solved and simplified using Mathematica V7 and the so- 
lutions are 

1 

aL = --48x- 
bL= l6 X 

(x 2 -lf ^ 

CR =" 



d.R 



Wx 
x 6 - 3x 2 + 2 
48x 



Finally the two parts of the Green's function are (introducing the x-dependence 
explicitly) 



4 2 x_ _ 1 i 

Mr,x) = -^- + ^- + ^-^-^ + 4-, l<r<x 
48x 16 r z 8 16a; 

, . . (x e -3x 2 + 2) (x 2 -l) 2 
Mr,x) = + *^<oc. 



(100) 



The solution for ^3 is then given formally by 

/r />oo 
jNL{x)(f>R{r,x)dx - J j NL (x)4> L (r,x)dx (101) 

This stream function is O(l) as r — > 00 so we claim that a secondary paradox of 
Stokes- Whitehead type has been avoided. To proceed further we need to better 
understand the full non-linear system. 



8 The non-linear equations and constraints on e 

The task now at hand is to develop a proper theoretical basis for the estimation 
of the parameter e. To this end we must write down the full Navicr Stokes 
equations. Converting to polar coordinates, the Navier-Stokes equations under 
consideration are 

r d(r,8) 
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Let us assume that R e is small enough that the flow remains symmetric about 
the horizontal axis. We may then write a Fourier decomposition 



^ = E sin(m0)(/> m (r) 



m— 1 



Then with Ck as above, 



VV = ^2 sin(m#)£ m m (r) 

m— 1 

oo 

VV = E sin(m0)O m (r) 



2 = 1 



XI sm(m9)C 2 m (t) m (r) 



m—l 



Re 

2r 



EE 

1=1 n=l 



sm(W)ncos(nO)— -£ n <j) n — lcos(W) sin(n#)0; — C 



oo oo r 

EE 

1=1 n=l 



dr 



rc[sin((Z - n)6) + sin((/ + n)6)]—C n (j) n 



dr 



- Z[sin((Z + n)6) - sin((Z - n)0)}^—C r , 



Doing the Fourier analysis gives us 

(i-<w) E(( m -^ £ 
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m — L t rm- 



or 



+E 



. Ci4>i + (l + m)(f)i +m ^Ci 



i=i 



dr 



(103) 

(104) 
(105) 



(106) 



(107) 



Thus far no approximations have been made, other than to assume the flow 
remains symmetric about the rc-axis. Now we approximate the model by con- 
sidering only the contributions of the terms with m = 1,2. We then have the 
simpler, but coupled and non- linear system 



2 Refd<f>i d 
C ^ 2= 2r\-drr Cl<t>1 ' <t>1 dr Cl4>1 



(108) 
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, 2 R e ( d(j) 2 „ ± ± d r ± n d(f>i „ ± , n± d 



Cifa = ^ ( -^rA& - </>i^£20 2 - 2^£ 2 02 + 202 — Ci<hJ (109) 
With our simple form of <pi that is a solution of the Laplace- Hclmholtz condition 

£ 2 l <f> 1 = e 2 £ l 4> 1 (110) 



we have already determined the solution of the Eq. ( 108 1- 02 is just ^3 as already 
given above. In order to generate maximal self-consistency of the linearization 
we must choose e in such a way as to minimize the mismatch between the right 



sides of Eq. (109) and Eq. ( 110 1 4 



Let us look at the overall scaling behaviour in terms of the Reynolds number. 
We let </> 2 = R e 4>2, so that 

The required matching is then of the form 

In a r2, K( d hr 1 1 d n 7 „ d <t>l n 7 , o7 d 



e £i<j)i ~ C x (j)i = — - — — C\4>\ - 017^^202 - 2— — £ 2 02 + 202tt£i^i 
2r \ or or Or or 

(H2) 

which strongly supports the scaling behaviour 

e 2 = Ai? 2 (113) 

as R e — > 0. That is, we have a theoretical basis for setting j3 = 1. A determina- 
tion of A requires the introduction of a suitable rigorous criteria for minimizing 
the mismatch between the linear and non-linear forms, and this is under inves- 
tigation. 

Note that the full Navier-Stokes equations, but limited to the first two an- 
gular modes, can certainly be written without further approximation as 

L\<\> x = ^A(r)A0! (114) 

for some unknown function A(r), and our method can now be properly under- 
stood as that of working with some "average" value of A(r), and noting that the 
resulting solution is free of a paradox. Note also that writing the right side of 
this linearized system as a multiple of L\<$>\ is not as arbitrary as it might seem, 
for we know that the right side of the full non-linear system vanishes identically 
when L\<$>\ = 0. That is, this Laplace-Helmholtz model has a proper theoretical 
justification, rather than merely being the basis of a convenient interpolation 
between the boundary conditions on the cylinder and at infinity. 

We could also consider generalizations where an improved ansatz for the 
form of X(r) is employed. The optimal average form and improved functional 
choices are under investigation. Only once this has been done would it make 
sense to consider further iteration. Some initial considerations suggest that as a 
function A(l) = and that A(r) might be asymptotic to a constant independent 
of R e as r — > oo, but further analysis is needed. 



4 This is of course one of many ways of proceeding, but our goal here is to provide a 
theoretical basis for the estimation of the so far free parameter e. 
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9 Summary 



We have constructed a global stream function satisfying the two-dimensional vis- 
cous incompressible steady Navier-Stokes equations in the limit R e — > 0. The 
stream function satisfies the correct boundary conditions on a cylinder and infin- 
ity. The "perturbation" to the linearized Navier-Stokes equations introduced 
to accomplish this is of higher order in R e than the inertia terms, and this 
modification is now properly understood as a linear representation of the full 
non-linear theory. The results suggest improved agreement with experimental 
data over those obtained by the MAE approach. Further work is needed on this 
approach, in particular on comparisons with newer data sets and determination 
of the remaining free parameter A from theoretical considerations. Although 
the approach of this paper initially started with rather ad hoc considerations, 
the method developed here is founded on a deeper consideration based on the 
approximation of the non-linear Navier-Stokes equations with an optimal linear 
approximation based on equations of Laplace-Hclmholtz type, rather than the 
biharmonic equation. The need to properly treat the non-linearities has also 
been illustrated by the observation that the presence of a paradox is unstable 
with respect to small changes in the system. 

9.1 Model summary 

The viscous and "paradox-free" stream function is given in non-iterated form 

by 

[ | 2Jr 1 ( g ) \l | 2Jr 1 (r g )' 



sin(0) (115) 



eK (e) ) r sK (e) 

and the parameter e = \/XR e . The drag coefficient with Tritton's conventions 
is: 

CD -R-lMef (U6) 

The choice of A from theoretical considerations has yet to be determined, but 
experimental drag data suggests A ~ 0.04452. The stream function satisfies the 
correct boundary conditions on the cylinder and at infinity for all e > 0, and 
also reduces to that for potential flow as e — > oo. Mathematica code for the 
stream function and velocity field are given in the Appendix. 

9.2 Credit 

I am grateful to Bin Zhou for his comments on the earlier (2006) version of this 
paper. 
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Appendix: Mathematica code for the flow field 

The following code may be useful. First the stream function: 

\[Psi][r_, \[Theta]_, \[Epsilon]_] := Sin [\ [Theta] ] (r - (1 + 

2 BesselK[l, \ [Epsilon] ] A [Epsilon] /BesselK [0 , \ [Epsilon] ] ) / 
r + 2 BesselK [1, r \ [Epsilon] ] A [Epsilon] /BesselK [0 , \ [Epsilon]]) 

The streamlines are then easily visualized: 

ContourPlot [\ [Psi] [Sqrt [x~2 + y~2] , ArcTan[x, y] , 0.01], {x, -10 , 
10}, {y, -5, 5}, AspectRatio -> 1/2, Contours -> 50, 
RegionFunction -> Function[{x, y}, x~2 + y~2 >= 1], 
Epilog -> Circle [{0, 0}, 1]] 



The Cartesian components of the velocity field are given by 
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Figure 2: Streamlines 



CartesianVelocity [x_ , y_, U_, \[Epsilon]_] := 

Module[{r = Sqrt[x~2 + y"2] , \ [Theta] = ArcTan[x, y] , Ur, U\[Theta]}, 
Ur = -((U*((-l + r~2)*\ [Epsilon] *BesselK[0, \[Epsilon]] - 

2*BesselK[l, \[Epsilon]] + 2*r*BesselK [1 , r*\ [Epsilon] ] ) * 
Cos [\ [Theta] ] ) / (r~2*\ [Epsilon] *BesselK [0 , \ [Epsilon] ] ) ) ; 
U\ [Theta] = 
U*(l + r"(-2) - (BesselK[0, 
r*\ [Epsilon]] - (2* 

BesselK[l, \ [Epsilon] ])/ (r~2*\ [Epsilon] ) + 
BesselK[2, r*\ [Epsilon] ]) /BesselK [0 , \ [Epsilon] ]) * 
Sin[\ [Theta]] ; 
Ur*{Cos [\ [Theta] ] , Sin [\ [Theta] ] > + 
U\ [Theta] *{-Sin [\ [Theta] ] , Cos [\ [Theta] ] >] 

The vector flow field is shown in Figure 3: 

VectorPlot [ 

If [x~2 + y"2 >= 1, 

CartesianVelocity [x, y, 1, 0.1], {0, 0}], {x, -10, 10}, {y, -5, 5}, 
RegionFunction -> Function[{x, y}, x~2 + y~2 >= 1], 
AspectRatio -> 1/2 , VectorScale -> Small, 
Epilog -> Circle [{0, 0}, 1]] 

This illustrates the satisfaction of the boundary conditions on the cylinder. 
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